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A  COMPUTATIONAL  MODEL  FOR  THE  SIMULATION  OF 
MILLIMETER-WAVE  PROPAGATION  THROUGH  THE  CLEAR  ATMOSPHERE 


Jerry  D.  Hopponen  and  Hans  J.  Liebe* 

Prediction  of  propagation  effects  (i.e.,  path 
attenuation,  phase  delay,  ray  bending,  and  medium  noise) 
over  the  1  to  300  GHz  frequency  range  through  the  clear, 
nonturbulent  atmosphere  is  accomplished  by  combining  a 
spectroscopic  data  base  with  a  computer  program  for  two- 
dimensional  ray  tracing.  Interactions  between  the  physical 
environment  and  electromagnetic  radiation  are  expressed  by  a 
complex  refractivity  N.  The  quantity  N  is  a  function  of 
frequency,  pressure,  humidity,  and  temperature. 
Spectroscopic  data  supporting  N  consist  of  more  than  450 
coefficients  describing  local  O2  and  H2O  absorption  lines 
complemented  by  continuum  spectra  for  dry  air  and  water 
vapor.  Height  profiles  (up  to  80  km)  of  N-spectra  are  the 
basis  for  calculating  propagation  effects  along  a  radio  path 
(ground-to-ground,  ground-to-aircraft ,  and  ground-to- 
satellite).  The  computer  model  assumes  a  symmetric,  spher¬ 
ically  sttatified  atmosphere  without  horizontal  N 
gradients.  Evaluation  of  path  integrals  for  radio  range, 
cumulative  attenuation,  and  noise  temperature  is 
accomplished  in  a  rapid  manner.  Various  simulated 
propagation  aspects  and  details  of  the  treatment  of  the 
noise  integrals  are  given. 

Key  words:  clear  atmosphere;  millimeter-wave  propagation;  path  attenuation 
and  delay;  radiances;  radio  path  modeling;  ray  bending 


1.  INTRODUCTION 

Performance  assessment  of  millimeter  wave  systems  that  transmit  or 
receive  through  the  Earth's  atmosphere  is  facilitated  by  considering  a  propa¬ 
gation  module  (climatological  influences)  and  a  nonpropagation  module  (hard¬ 
ware,  software,  etc.).  For  reliable  system  design  as  well  as  for  adaptive 
compensations  of  atmospheric  propagation  effects,  it  is  imperative  that  a 
model  for  simulating  atmospheric  transfer  characteristics  be  available.  The 
formulation  presented  here  combines  a  broadband  model  of  the  complex  refrac¬ 
tive  index  of  clear  air  with  a  ray  tracing  program  in  order  to  calculate  radio 
path  properties  over  the  extended  EHF  (1  to  300  GHz)  frequency  range. 


*Jerry  D.~  Hopponen  is  with  Lockheed  Missiles  and  Space  Company,  Sunnyvale, 
CA  94086  (formerly  with  NT IA/ ITS,  Boulder,  CO);  Hans  J.  Liebe  is  with  the 
Institute  for  Telecommunication  Sciences,  National  Telecommunications  and 
Information  Administration ,  U.S.  Department  of  Commerce,  325  Broadway, 
Boulder,  CO  80303-3328. 


The  clear,  nonturbulent  air  mass  up  to  an  altitude  of  80  km  acts  as  a 
unique  filter  of  the  EHF  band  with  attenuation,  delay,  and  noise  char¬ 
acteristics  not  found  at  frequencies  be^w  20  GHz.  Absorption  of  radio  waves 
by  molecular  oxygen  and  water  vapor  results  in  frequency-dependent 
attenuation,  delay,  and  atmospheric  emission;  all  of  which,  in  return,  offer 
remote-sensing  opportunities  on  the  atmospheric  state.  Propagation  effects 
depend  on  the  climatology  assumed  for  the  Earth-to-space  path  and  may  show 
considerable  diurnal  and  seasonal  variation.  Therefore,  although  standard 
reference  values  may  be  available  for  the  phenomena  of  interest,  the  careful 
assessment  of  propagation  effects  requires  a  computer  model  capable  of  making 
accurate  predictions  using  any  set  of  meteorological  data. 

The  interaction  of  propagating  radio  waves  with  the  atmosphere  is 
described  by  a  complex  refracti vity,  N.  In  the  EHF  band,  the  frequency  depen¬ 
dent  components  of  N  are  primarily  due  to  the  absorption  spectra  of  molecular 
oxygen  and  water  vapor.  Theoretical  (Rosenkranz,  1975;  Smith,  1981)  and 
experimental  (Liebe  et  al.,  1977;  Liebe,  1983;  and  references  in  Liebe,  1985) 
efforts  have  produced  analytical  expressions  that  are  used  in  conjunction  with 
a  spectroscopic  data  base  for  the  calculation  of  N  as  a  function  of  frequency 
and  meteorological  variables.  Characteristics  of  EHF  propagation  are 
simulated  by  generating  a  profile  of  N  from  a  height  distribution  of  clear  air 
data  and  employing  the  N-profile  in  a  ray-tracing  methodology.  Some  details 
of  the  physical  models  and  algorithms  incorporated  in  the  complete  simulation 
program  are  discussed  in  order  to  set  up  a  computer  program  in  optimum 
manner . 

2.  ATMOSPHERIC  INTERACTION  WITH  RADIATION 

The  first  step  in  a  simulation  of  electromagnetic  wave  propagation 
through  the  first  80  km  of  the  Earth's  atmosphere  is  a  specification  of  the 
state  of  the  medium.  For  a  clear,  nonturbulent  atmosphere  it  is  sufficient  to 
give  height  profiles  of  total  pressure,  P;  ambient  temperature,  T;  and 
absolute  humidity,  v.  These  quantities  are  described  by  data  specified  (as  in 
reference  atmospheres  or  radiosonde  data)  at  discrete  altitudes  all  expressed 
in  standard  SI  units  (T  in  degrees  Kelvin,  P  in  Pascals,  v  in  grams  per  cubic 
meter,  and  height  in  meters).  The  requisite  height  transformation  is  given  by 
List  (1958)  and  the  spectroscopic  relations  are  found  in  Liebe  (1985). 


Properties  of  the  propagation  medium  are  expressed  by  a  complex  refract- 
ivity  N( P,T,v) .which  consists  of  three  components 

N  =  N0  +  D(f )  +  j  N" (f )  ppm  (1) 

where  N0  is  the  frequency-independent  refracti vity ,  and  D(f)  and  N"(f)  are 
frequency  dependent  dispersion  and  absorption  terms  arising  from  spectra 
exhibited  by  molecular  oxygen  and  water  vapor.  The  N-formul ation  is  founded 
on  up-to-date  theoretical  results,  incorporates  empirical  correction  terms, 
and  is  amenable  to  rapid  execution  on  a  main-frame  computer. 

Amplitude  and  phase  responses  of  a  plane  radio  wave  generally  are 
described  by  the  propagation  constant 

T  =  -(a/2)  +  j (20958f  +  0)  km-1  (2) 

where  the  specific  power  attenuation  is 

a  =  0.04192fN"  Np/km,  (3a) 

or,  in  units  of  decibel  per  km, 

a*  -  4.343a  dB/km;  (3b) 

the  specific  phase  delay  is 

B  =  0.02096f(NQ  +  D)  radians/km  (4) 

or,  expressed  in  excess  propagation  delay, 

t  -  (8/2nf)103  =  3. 336(Nq+  D)  ps/km;  (5) 

and  the  frequency,  f  is  in  units  of  gigahertz  throughout.  In  a  nonhomogeneous 
medium  the  variable  N  is  integrated  along  a  ray  path,  as  discussed  below.  The 
complex  refracti vity ,  N,  expressed  in  terms  of  measurable  quantities,  provides 
a  macroscopic  measure  of  the  interaction  between  radiation  and  the  absorbing 
molecules  in  moist  air. 

Prediction  of  the  propagation  effects  along  short  horizontal  paths  may  be 
accomplished  via  expressions  (2)  through  (5).  For  the  case  of  a  curved  path 


refracted  through  a  nonhomogeneous  clear  air  mass  (e.g.,  ground-to-satel 1 ite 
links,  radar  tracking,  remote  sensing),  meteorological  data  profiles  are  used 
as  input  and  then  converted  to  corresponding  profiles  of  N  values  for  use  in 
ray-tracing.  It  is  advisable  to  store  the  N  profiles  on  a  computer  disk  file 
since  ray  tracing  may  be  performed  by  several  program  executions.  This 
approach  affords  some  economy  since  it  permits  a  separate  main  program  for 
calculating  N.  Reducing  the  compiled  size  of  the  ray-tracing  program  embraces 
a  philosophy  of  modular  programming.  It  may  happen  that  the  atmospheric  data 
are  specified  at  altitudes  that  are  relatively  far  apart,  so  that  the  derived 
N  profile  is  of  questionable  value  in  ray  tracing.  Thus,  interpolation  on 
either  the  basic  atmospheric  data  or  on  the  N  profile  is  required.  Although 
working  with  the  atmospheric  data  that,  in  a  sense,  count  the  absorber 
molecules  of  the  path,  there  are  advantages  to  having  more  height  resolution 
directly  in  the  N  profile.  Each  addition  to  the  atmospheric  data  profile 
requires  three  interpolations  (for  T,  P,  and  v),  which  would  be  followed  by  a 
calculation  of  N.  Such  an  approach  is  not  only  time  consuming  but  also 
presumes  that  the  ray-tracing  program  has  the  code  to  calculate  N,  contrary  to 
the  philosophy  advocated  earlier.  Hence  the  necessary  ''closeness"  of  interpo¬ 
lated  N  values  is  determined  by  a  convergence  criterion  in  the  adaptive  numer¬ 
ical  integration  employed  to  evaluate  the  ray-tracing  integrals.  The  number 
of  altitudes  needed  is  generally  not  known  prior  to  program  execution  and 
depends  on  the  given  atmospheric  data  and  the  ray  elevation  angle. 

The  N(h)  structure  over  the  first  3  km  above  surface  level  hQ  is  critical 
for  a  successful  prediction  of  refraction,  affecting  rays  starting  at  hQ  under 
low  ( i .  e . ,  close  to  the  horizontal)  elevation  angles,  0Q.  Values  of  N ( h )  that 
are  missing  in  the  input  data  are  generated  by  interpolation.  The  inter¬ 
polation  for  the  real  part  of  (1)  neglects  the  dispersive  term  0(f) 

(generally,  D  <  0.01No)  and  assumes  an  exponential  decay  of  moist  air  refrac- 
tivity  with  height  in  the  form 

NQ(h)  =  No(ho)exp^ho  '  h)/hs3  ppm-  (&a) 

Surface  refractivity  and  scale  height  h$  are  interrelated  on  the  average 
via  (8ean  and  Thayer,  1959) 

h^  =  1/ In  Ll / LT _a ]  km  (6b) 
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where 


a  =  (7.32/Nj)  exp  (0. 00558N^) .  (6c) 

Typically,  hs  varies  between  8.0  and  4.5  km  when  changes  from  280  to 
450  ppm  depending  on  surface  humidity. 

Linear  interpolation  for  values  of  N"(f)  can  be  obtained  at  heights  other 
than  those  of  the  original  data  profile.  Ray-tracing  calculations  employing 
exponential  interpolation  yield  values  that  differ  at  select  frequencies  that 
differ  by  less  than  3  dB  from  those  calculated  with  linear  interpolation  when 
the  lower  20  km  of  atmospheric  data  are  given  at  1-km  altitude  intervals. 
However,  it  should  be  noted  that  when  the  same  data  set  is  reduced  to  model 
5  km  intervals,  the  disparity  between  the  exponential  and  linear  ray  trace 
results  grows  to  over  30  dB.  Specific  results  are  shown  below  in  Table  1. 

The  scale  height  for  absorption  varies  with  humidity  distribution  and 
frequency  (e.g.,  between  11  and  1.9  km;  see  Liebe,  1985:  TABLE  4). 

By  virtue  of  the  interpolation  methods  for  the  real  and  imaginary  parts 
of  N,  the  spherically  symmetric  shells  defined  by  the  N  profile  can  be 
subdivided  into  thinner  shells,  thus  more  nearly  approximating  a  continuum  of 
N  values. 


TABLE  1. 

Interpolation  Comparison  for  a  Tangential  Path  (0q  =  0°)  Through  the  U.S. 


Standard  Atmosphere  (h 

«  0  to  80  km) 

FULL  DATA 

SET 

SPARCE  DATA 

SET 

Ah  »  1  km 

Ah  -  5  km 

Frequency 

Attenuation  A (d B) 

Attenuation  A(dB) 

(GHz) 

Linear 

Exponential 

Linear 

Exponential 

20 

34.66 

34.24 

40.63 

32.51 

40 

36.57 

36.16 

43.33 

36.75 

60 

5749.7 

5744.8 

5803.1 

5722.5 

95 

117.82 

115.87 

146.12 

112.31 

120 

491.99 

490.04 

537.89 

508.07 

5 
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3.  RAY  PATH  CHARACTERISTICS 


3 . 1  Ray  Tracing 

Wave  propagation  over  a  radio  path  can  be  described  when  increments  dx  of 
the  total  path  length  i,  both  in  km,  are  assumed  in  the  x  direction. 
Simplifications  are  introduced  that  allow  the  application  of  numerical  methods 
to  three  types  of  radio  paths: 


Path  Type 

Height  Range 

Path  Differential 

h 

dx 

hori zontal 

i(<  50  km) 

verti cal 

hQ  to  H 

Ah 

si  ant 

hQ  to  H 

(Ah)s 

The  nonhomogeneous  (vertical,  slant)  atmosphere  is  assumed  to  be  spherically 
stratified  in  concentric  layers  of  height  intervals  Ah,  for  which  values  of  N 
can  be  specified. 

Ray  tracing  is  used  for  slant  paths  to  determine  the  path  extension 
factor  defined  by 


s  *  dx/dh  (7) 

A  "ray"  may  be  regarded  as  a  curve  generated  by  successive  normals  to  surfaces 
of  constant  phase  of  a  propagating  wave.  For  a  nonhomogeneous  atmosphere  the 
phase  velocity  may  vary  from  point  to  point,  resulting  in  ray  bending.  A  ray 
is  assumed  to  start  from  the  initial  level  hQ  with  an  elevation 
angle  0Q  (measured  from  the  horizontal)  and  proceeds  through  the  atmosphere 
gaining  the  height  interval  Ah  =  (h  -  h  ),  while  being  subjected  to  refractive 
bending  when  0Q  <  90°.  During  stable  climatic  conditions  (i.e.,  a  homogene¬ 
ously  stratified  atmosphere),  the  path  extension  of  a  slant  path  down  to 
about  0Q  i  10°  increases  according  to  the  secant  law 
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Zenith  (0  =  90°)  path  behavior  is  said  to  be  caused  by  "one"  air  mass. 

The  effects  of  both  refraction  and  Earth's  curvature  determine  for  low- 
angle  (0  <  10°)  cases  the  path  increments  dx  =  (Ah)s,  which  are  given  by 
(Blake,  1968;  Levine  et  al  . ,  1973) 


- 1  fe)  m-  *•] 


where  =  [NQ  +  D(f)Jh  ,  Nj  =  [N0  +  0(f) ]h ,  and  r^  =  6357  km  is  a  standard 
(45°  N)  Earth  radius.  Equation  (7b)  is  the  spherical  form  of  Snell's  law  and 
is  strictly  valid  only  when  N  is  purely  real.  In  the  case  at  hand,  values  for 


N"  are  generally  small  compared  with  N- 


For  very  low  elevation 


angles  (0  <  1°)  and  values  of  h  close  to  hQ,  many  fine  steps  of  Ah  are 
required  and  (7b)  becomes  numerically  unstable.  Approximations  of  the  two 
terms  containing  j  and  r^  are  made  by  expanding  both  as  power  series  and 
retaining  terms  of  degree  less  than  three  (e.g.,  Blake,  1968).  At  this  point, 
increments  dx  =  (Ah)s  become  very  sensitive  to  the  height  distributions  of 
meteorological  variables.  For  example,  the  tangential  (0Q  =  0°)  air  mass 
through  the  U.S.  Standard  Atmosphere  for  dry  air  is  38  (35)  times  the  zenith 
value,  while  for  water  vapor  the  tangential  mass  can  vary  betwen  70  (64)  and 
180  times  the  zenith  case.  Values  in  parenthesis  give  the  air  mass  for 

=0  (that  is,  no  refractive  bending). 

The  path  differential,  dx,  is  used  in  the  formulation  of  several 
integrals  of  interest.  These  integrals  resist  attempts  at  closed-form  solu¬ 
tion,  owing  to  the  height  dependence  in  dx,  and  so  must  be  evaluated  by  numer¬ 
ical  methods.  The  adaptive  numerical  integration  technique  by  Romberg  is 
found  to  provide  accurate  and  rapid  results  because  of  the  way  in  which  the 
input  data  is  treated  and  also  because  of  the  nature  of  the  integrands  (Bauer, 
1961).  The  Romberg  method  not  only  inserts  additional  points  into  a  given 
height  interval  but  also  shifts  to  higher  order  schemes  [e.g.  trapezoidal  rule 
(linear)  to  Simpson's  rule  (quadratic)  and  so  on],  thus  accelerating  the 
convergence  of  the  algorithm.  The  adaptive  feature  is  especially  desirable 
since  the  number  of  points  required  to  evaluate  an  integral  across  an  atmos¬ 
pheric  layer  may  vary  as  a  function  of  layer  "thickness"  (i.e.,  the  difference 
between  the  initial  and  terminal  altitudes  defining  the  layer)  and  the  ray 
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in  the  lowest  layer,  so  this  layer  is  always  subdivided  into  groups  of  10 
intervals  of  thicknesses:  0.5,  5.0,  50,  .  .  .  meters  up  from  the  surface  to 
the  first  data  point  above  the  surface.  These  artificial  layers  are  taken  as 
the  intervals  of  integration  in  the  lowest  part  of  the  atmosphere. 

3.2  Path  Integrals 

Radio  range  (real  part  of  N) 


R  -  /  [N0  (h)  ♦  D(f,h)]dx  mm,  (8) 

ho 

particularly  for  low  elevation  angle  cases,  0Q  <  5°,  is  the  first  path 
integral  to  be  evaluated  by  means  of  the  ray-tracing  routine.  The  dispersion 
term,  D(f),  adds  a  f requency-dependent  component.  Rays  propagating  at 
different  frequencies  will  separate  in  space  (and  time).  A  radio  range 
difference,  AR,  can  be  calculated  by  tracing  rays  with  and  without  dispersion, 
D (f ) .  The  example  in  Figure  1  for  a  tangential  path  (0Q  -  0)  reveals  differ¬ 
ences  on  the  order  AR  -  0.3  and  -0.4  m  at  the  frequencies  57  and  63  GHz, 
respectively.  In  theory,  the  two  rays  are  separated  by  Ah  -  0.7  m  when 
exiting  the  atmosphere  at  H  =  80  km;  in  practice,  path  attenuation  A  [see  (9) 
below]  on  the  order  to  6000  dB  disallows  this  case. 

Digital  transmission  at  high  data  rates  (>0.2  Gbit/s),  however,  will  be 
impaired  at  carrier  frequencies  where  the  atmosphere  displays  gradients  dD/df 
across  the  channel  bandwidth  (>0.5  GHz).  An  example  of  pulse  propagation 
(1  ns,  trapezoidal  shape)  over  a  horizontal  (homogeneous)  path  is  shown  in 
Figure  2.  Shape  distortion  and  pulse  widening  increase  with  the  propagation 
distance  1.  The  time  domain  analysis  (K.  Allen  and  H.  Liebe;  private  communi¬ 
cation,  1985)  has  not  been  extended,  so  far,  to  address  pulse  transmission 
affected  by  R  (and  A) . 

Total  attenuation  along  the  ray  path  (imaginary  part  of  N)  is  given  by 


J  o' 


(f ,h)dx 


where  a*(h)  is  specific  attenuation  in  dB/km  (3b)  as  a  function  of  height. 
Attenuation  A,  quantifies  the  amount  of  energy  extracted  from  a  plane  wave  for 
a  one-way  path  through  the  atmosphere.  Path  attenuation  may  show  considerable 
variations  with  season,  as  depicted  by  two  examples  in  Figure  3 
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Tangential  Path  Dispersion 
(Altitude  Dependence) 
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Figure  1. 


Radir  range  difference  AR  -  R ( N  +  D)  -  R(NQ)  due  to 
atmospheric  phase  dispersion  D(f,h)  for  the  frequency  pair 
57/63  GHz  propagating  at  an  elevation  angle  of  0°  (tangential 
path)  from  sea  level  (h=  0  km)  through  a  midlatitude  summer 
atmosphere  (Valley,  1976).  The  refractive  radio  range  (H  =  80 
km)  is  R(N0)  =  103.74  m;  path  attenuation  A  is  on  the  order  of 
6000  dB  at  both  frequencies  (Liebe,  1983). 


Pulse  Distortion  at  Sea  Level 
(Path  Length  Dependence) 


PULSE  WIDTH,  ns 

Figure  2.  Normalized  (amax  =  1)  trapezoidal  1  ns-pulse  at  the  carrier 
frequency  fc  =  55  GHz  propagating  over  a  terrestrial  path  at 
sea  level  various  distances,  l  =  0  to  20  km. 

The  following  meteorological  conditions  are  assumed: 

P  =  1013  mb, 

T  =  1 5°  C , 

v  =  11.5  g/m3  (90%RH) , 

leading  at  fc  =  55  GHz  to  the  specific  values 
a*  =  4.29  dB/km  and 

B  =  0.916  radians/km. 
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for  a  dry  midlatitude  winter  atmosphere  and  a  moist  summer  atmosphere.  The 
fine  structure  resulting  from  the  complex  of  O2  lines  centered  around  60  GHz 
is  not  shown  (see  Liebe,  1983),  but  a  crossover  of  the  two  curves  is 
indicated,  which  results  from  the  complicated  temperature  dependence  in  the 
spectral  absorption  expressions. 

Radiance  Tg  indicates  noise  emission  so  that  the  atmosphere  remains  in 
thermal  equilibrium  and  is  expressed  as  equivalent  blackbody  temperature  of 
radiation  by  the  Ray  1 ei gh-Jeans  approximation  to  Planck's  law  (Chandrasekar , 
1950)  in  the  form 

Tb  =  f T(x)  $(x)dx  K,  (10) 


where 


*(x)  =  a(x)t(x‘ ,x) 


(10a) 


is  a  weighting  function  on  T(x),  the  local  ambient  temperature  along  the 
path;  x'  and  x  are  two  adjacent  points;  and  t  is  a  transmittance  to  be  defined 
below.  Two  cases  can  be  made.  In  the  first  one,  radiation  upwelling  from  a 
starting  height  hQ  to  the  outer  limit,  H  of  the  absorbing  atmosphere  (and 
possibly,  continuing  to  a  nadir-viewing  satellite)  leads  to 


Tj  -  /  *uT(x)dx 

h 

0 


(11a) 


(surface  emission  is  not  considered).  Secondly,  radiation  downwelling  to 
height  hQ  (not  necessarily  the  Earth’s  surface  height)  is  described  by 


T*  -  /  V(x)dx 


(lib) 


(cosmic  oackground  radiation  of  *2.7  K  is  not  considered).  The  correspondi  ng 
transmi ttances  in  $  are  defined  by 


x  (h,H)  =  exp[  -J  o(x)dx] 
u  If 


(12a) 


Wh)  =  exp^  ".J a(x)dx] . 
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Figure  3. 


Seasonal  variation  of  total,  one-way  attenuation  A  through  a 
dry  winter  and  moist  summer  midlatitude  atmospheres  (Valley, 
1976)  along  a  vertical  (zenith)  Earth-to-space  path. 


They  give  the  power  loss  from  a  point  of  height  h  to  the  reference  heights,  H 
or  hQ.  In  both  integrals,  dx  is  a  differential  of  curved  ray  path  (7b).  The 
factor  aT  is,  by  the  Raleigh-Jeans  approximation,  proportional  to  emission 
from  a  volume  element  of  the  absorbing  air. 

Figure  4  shows  atmospheric  noise  due  to  vertical  radiation  upwelling  to  a 
satellite  through  two  atmospheres  representing  extreme  environmental 
conditions.  The  moist,  tropical  air  mass  generally  yields  higher  noise 
temperatures  than  the  dry,  subarctic  winter  atmosphere.  The  actual  behavior 
is  dependent  on  the  selected  frequency  as  well  as  the  physical  conditions 
between  0  (hQ)  and  80  (H)  kilometers. 

Figure  5  compares  atmospheric  contributions  to  upwelling  and  downwelling 
radiances  through  a  vertical  path  in  a  midlatitude  winter  atmosphere.  Except 
for  regions  of  strong  absorption,  where  the  ground-based  sensor  receives 
radiation  only  from  the  immediate  vicinity,  the  curves  are  nearly  coincident 
across  the  EHF  band. 

Weighting  functions  $u  peaking  at  58.82  GHz  in  the  stratosphere  are 
exhibited  in  Figure  6.  The  elevation  angle  0Q  at  hQ  =  0  varies  the  height 
from  where  maximum  radiance  originates  between  h  -  18  (Qq  «  90°)  and 
h  -  25  km  (0Q  -  0°).  The  radiance  that  can  be  detected  under  these  conditions 
at  H  =  80  km  (and  beyond)  is  displayed  in  Figure  7.  No  radiance  is  received 
from  the  ground  since  the  zenith  attenuation  A(0Q  =  90°)  at  this  frequency  is 
140  dB.  The  temperature  height  profile  of  the  particular  model  atmosphere 
accounts  for  the  radiance  structure  as  suggested  by  Figure  6.  At  this 
frequency  there  is  less  than  2  deg  K  variation,  so  the  Earth  appears  from 
outer  space  as  a  nearly  uniform  noise  source. 

4.  NOISE  TEMPERATURE  CALCULATION 

The  ray  tracing  methodology  that  is  employed  to  compute  the  noise  temper¬ 
atures  above  begins  with  the  specification  of  atmospheric  data  at  heights 

ho  <  hi  <  h2  .  .  . 

The  noise  temperature  integrals  can  then  be  written  as 
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NOISE  TEMPERATURE 


Zenith  Path  Radiance 
(Climatological  Extremes) 


EREQUENCY  f,  GHz 


Figure  4.  Upwelling  atmospheric  noise  temperature  (radiance)  Tn  from 

tropical  and  subarcti c  model  atmospheres  (Valley,  19/6)  repre¬ 
senting  environmental  extreme  conditions  for  a  vertical  Earth- 
to-space  path. 


NOISE  TEMPERATURE 


Zenith  Path  Radiance 
(Upwelling/Downwelling) 


Figure  5.  Atmospheric  noise  temperatures  Tg  for  upwelling  (no  surface 
contribution)  and  downwelling  (no  2.7  deg  K  cosmic  background 
contribution)  radiation  along  a  vertical  path  through  a  mid¬ 
latitude  winter  atmosphere  (Valley,  1976). 
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Upwelling  Radiance  at  58.8  GHz 
(Elevation  Angie  Dependence) 
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Figure  7. 


ELEVATION  ANGLE  0O  ,  deg 


Upwelling  atmospheric  noise  temperature  TB  at  58.82  GHz  leaving 
the  atmosphere  at  H  »  80  km  for  rays  that  originate  at  h  -  0  km 
under  various  angles  0  *  0  to  90s  (model  atmosphere,  see 
Figure  6). 
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*  m  /*h,.  i 

t’b  -  L  Tu(hitl,H)  /  1*\(h.h1ti)oTdx 


(13a) 


T8  '  hMP'1  Td(h1  ,h)o,Td> 


(13b) 


In  the  expression  for  the  upwelling  temperature,  the  index  m  is  such  that 
hm+l  =  ^ >  wh^e  f°r  the  downwelling  case,  n  is  such  that  hp+j  is  a  suitable 
replacement  for  the  infinite  limit  of  integration  which  appears  in  (lib). 

This  is  not  only  convenient  but  also  physically  meaningful  since  above 
approximately  80  km,  the  absorption  is  nearly  nonexistent  (Liebe,  1983)  and 
the  assumption  of  thermodynamic  equilibrium  begins  to  break  down  (Smith, 

1982) . 

Although  the  summands  in  the  equations  for  upwelling  and  downwelling 
noise  temperature  can  be  calculated  and  summed  by  the  ray  tracing  program,  a 
pair  of  recursion  relations  that  actually  yield  noise  temperatures  at  not  only 
hQ  and  H,  but  also  at  all  other  heights  h^  ,  are  preferable  to  direct 
summation.  For  noise  temperature  due  to  radiation  flowing  along  a  ray  up  to 
height  hk+1,  it  follows  that 

rh  k+1 

TB<hk+l>  -J  Vh’hk+l>aTdx  +  T>k)Vhk’hk+l)  K-  (14a> 


Similarly,  for  downwelling  radiation 


w 


/*k+ 1 

- J  Tn(hu» 


d(hk,h)aTdx  +  TB(hk+1)xd(hk,hk  +  1) 


(14b) 


Heuri sti cal ly ,  the  recursion  relations  show  that  the  noise  temperature  at  one 
side  of  a  layer  is  the  sum  of  the  emission  of  that  layer  plus  the  noise 
temperature  at  the  opposite  side  diminished  by  a  transmission  factor. 

In  view  of  these  equations,  it  is  sufficient,  in  the  course  of  a  surface- 
to-space  ray  trace,  to  compute  and  store  the  quantities 


Tu(hk*hk+1^ 


k  +  1  /*^k  +  l 

J  Tu(h,hk+1)aTdx,  and  J  xd(hk, 


h)aTdx 


and  afterwards  to  combine  them  as  indicated  above. 


If  the  temperature  T 1  at  the  altitude  is  perturbed  pwhen,  for  example, 
it  is  required  to  obtain  the  Jacobian  matrix  of  partial  derivatives  in  the 
Backus-Gi 1 bert  inversion  method  (Westwater  and  Strand,  1974)],  only  those 
integrals  with  h^  as  an  upper  or  lower  limit  need  be  recomputed,  and  the 
recursion  relation  then  applied. 

The  calculation  of  total  attenuation  across  a  layer  is  direct.  The  noise 
temperature  integrals  across  a  layer  may,  in  principle,  be  evaluated  by  the 
adaptive  Romberg  subroutine  (Bauer,  1961)  which  will  require  values  of  the 
i ntegrands 

Vh,hk+l)a(h)T(h)  dnd  Td(hk  >h)<*(h)T(h) 

at  a  number  of  intermediate  heights  between  h^  and  h|<  +  ^.  The  exact  number  of 
intermediate  heights  is  determined  as  a  function  of  the  convergence  criterion 
in  the  adaptive  Romberg  integration  routine. 

Since  the  evaluation  of  each  integrand  involves  the  computation  of 
attenuation  across  a  sublayer  of  the  original  layer,  such  a  computation  is,  in 
a  sense,  redundant.  Furthermore,  computing  time  increases  (and  so  does 
expense)  as  the  number  of  calls  to  the  integration  routine  increase.  Finally, 
the  exponential  factors  in  the  noise  integrands  are  nearly  the  same,  except 
that  one  goes  from  the  top  of  the  layer  downward,  while  the  other  goes  from 
the  bottom  up,  thereby  introducing  yet  another  apparent  redundancy. 

Therefore,  in  the  interest  of  more  efficient  computation,  several  numerical 
methods  have  been  devised. 

It  is  first  noted  that  if  the  temperature  T  is  constant  between  two 
successive  heights,  a  and  b  (which  may  lie  between  hk  and  h^)  then 

fb 

xu(h,b)aTdx  =  T|_l  -  x(a,b)]  =  J  x^ (a ,h )aTdx .  (15) 

d 

This  equation  provides  a  simple  means  of  evaluating  the  noise  integrals  across 
an  isothermal  layer,  as  is  often  postulated  or  measured  in  the  upper  air 
mass.  When  the  temperature  between  h^  and  h|<  +  i  is  not  constant,  the  layer  may 
be  subdivided  into  small  subintervals  such  that  the  temperature  over  each  is 


nearly  constant.  In  that  case,  the  noise  integrals  across  the  layer  are 
approximated  as 


(h,hk+1)aTdx  = 


E  ‘  Wam>]  >16a> 


and 


,h)aTdx  = 


E  WW1  -  VWl'J-  t16b» 


Here  the  points  aQ  =  h^  <  a^  <  a2  .  .  .  <aL+i  =  hk=1  define  the  subintervals 
of  the  layer  between  h^  and  hk=^  and  the  number  of  points  a-j  is  chosen  so  that 
the  temperature  is  nearly  constant  between  a^  and  a^+^  for  all 
i  -  0,  I,  ...  L.  Each  summand  is  readily  computed  once  the  intergrals 

i  +1 

a  dx,  i  -  0,  1,  2,  ...  L 

are  available.  Although  these  may  be  calculated  by  calling  the  integration 
subroutine,  a  more  efficient  approach  can  be  devised. 

Between  heights  hk  and  hk+1  a  curved  ray  path  has  the  length 


f^k 


where  dx  =  s ( Ah )  is  the  differential  path  length  (7b).  If  the  layer  is  parti¬ 
tioned  into  L  +  1  subintervals,  as  above,  then  a  differential  of  path  length 
across  the  interval  from  ai  to  ai+^  is  approximately 

dP  -  -lT r,  (18) 

and  if  there  is  no  ray  bending  then  each  subsegment  has  length  equal  to  dP. 

When  the  elevation  angle  0  is  assumed  constant  and  the  attenuation  varies 
slowly  between  h k  and  hk+j  then  dx  =  dh/sinQ  and 


a(h)dx  = 


a(ai+1)  ♦  a(a.)  a.+1 


-  a 


a(ai+1)  +  a(ai)  r  U\ 


sin  0 


/ 

a . 

l 


Let  denote  the  average  value  of  a(h)  between  a^  and  a^,  so  that,  if 
constant,  it  follows  from  the  definition  of  dP  that 

i  +1 

adx  =  a^dP 


dx . 


0  is 


and  if  0  is  not  constant  then  this  is  an  approximation.  Consequently,  the 
difference 


A+i 

= J  adx  -  I]  a, 
h,  i =0  1 


may  not  equal  zero.  To  retain  the  simplicity  of  the  summation,  a  correction 
term  for  the  attenuation  values  is  defined  as 


Then 


as  desired. 
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A 

( L+l )  dP  ' 


£0  (a,  *  «)dP  -  A  -  E-(uTTdrdP  ■ 0 
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5.  CONCLUSIONS 


Computer  simulation  of  propagation  effects  is  essential  to  the  design  of 
systems  operating  through  the  atmosphere.  Molecular  absorption  varies  across 
the  1  to  300  GHz  oand  and  poses  constraints  for  transmission  and  shielding 
applications.  Transfer  properties  are  strongly  dependent  upon  radio  path 
length  and  elevation  angle  and  also  on  meteorological  conditions.  Methods 
have  been  presented  to  model  these  properties  based  on  physical  principles, 
algorithm  stability  and  efficiency,  and  computing  philosophy.  The  noise 
temperature  algorithms  not  only  afford  a  means  of  obtaining  noise  temperatures 
at  many  points  along  a  ray  (only  at  the  expense  of  the  time  required  to  print 
the  additional  information)  but  also  provide  an  elementary  method  for 
obtaining  the  elements  of  the  Jacobian  matrix  used  in  temperature  inversion. 
The  major  limitations  of  the  model  result  from  the  great  variability  of  the 
atmosphere  and  the  sporadic  occurrence  of  clouds  and  rain,  which  were  not 
addressed , 
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